Granular flow down an inclined plane: Bagnold scaling and rheology 
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We have performed a systematic, large-scale simulation study of granular media in two- and three- 
dimensions, investigating the rheology of cohesionless granular particles in inclined plane geometries, 
i.e., chute flows. We find that over a wide range of parameter space of interaction coefficients and 
inclination angles, a steady state flow regime exists in which the energy input from gravity balances 
that dissipated from friction and inelastic collisions. In this regime, the bulk packing fraction (away 
from the top free surface and the bottom plate boundary) remains constant as a function of depth 
z, of the pile. The velocity profile in the direction of flow Vx{z) scales with height of the pile H, 
according to Vx{z) oc H" , with a = 1.52 ± 0.05. However, the behavior of the normal stresses 
indicates that existing simple theories of granular flow do not capture all of the features evidenced 
in the simulations. 
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I. INTRODUCTION 



It is tempting to regard the behavior of granular ma- 
terials as being a problem in engineering or applied sci- 
ence, inasmuch as the fundamental laws governing their 
constituent particles are well known. Being comprised 
of macroscopically large grains, granular materials obey 
classical mechanics, although the existence of friction and 
inelastic collisions complicates matters. However, while 
it is true that the collision of two grains is analytically 
tractable, an aggregate of such grains is a many-body 
system, whose macroscopic behavior cannot be simply 
related to the laws controlling individual constituents. 

For this reason, a continuum treatment is often 
adopted, in which the variables are averaged properties 
whose governing equations are derivable, in principle, 
from the known microscopic laws. Among these aver- 
aged variables are the density p and the stresses aaf3, 
which obey the Cauchy equations that enforce momen- 
tum conservation (or force balance if there are no accel- 
erations). However, this set of equations is insufficient to 
solve for the stresses, since there are too few equations: 
in D dimensions, there are D(D + 1)/2 independent com- 
ponents of (Tq/3 (which is a symmetric tensor), but only 
D equations of momentum conservation. Therefore, the 
Cauchy equations must be augmented by additional con- 
stitutive relations, possibly history-dependent, which tell 
how the material in question responds to the application 
of a force. It is in these constitutive relations that the 
specifics of the material in question come into play. In 
the case of steady state flow, which we will consider in 
this paper, constitutive equations would relate the strain- 
rate, 7q,/3, to the stress. 

In 1954, Bagnold 0] proposed that in inertial granular 
flow, the shear stress is proportional to the square of the 



strain-rate: 

a oc 7^. (1) 

His argument, applied to the case of bulk granular flow, 
is predicated on a constant density profile. In practice, 
the presence of significant finite-size or wall effects often 
obscures Bagnold scaling. In this study, we report a set 
of numerical simulations of bulk granular fiow down an 
inclined plane, so-called "chute fiow" , in two and three 
dimensions. The geometry is simple: a layer of bulk gran- 
ular material is placed on a fiat plane of area A (or line 
of length L in 2D) on which grains have been glued, so as 
to form a rough base. The thickness of the layer is mea- 
sured in terms of the pile height parameter H = N(P /A 
(or Nd/L in 2D), where N and d are the total number 
of particles and their diameter, respectively. The plane 
is inclined at an angle 9 and the flow is observed. The 
parameters controlling the fiow are the macroscopic vari- 
ables 6 and H, as well as the microscopic variables de- 
termining the nature of interaction between two grains, 
such as grain friction ^ and coefficient of restitution e. 

In Ref. 1^ , we provided a summary of our simulations 
in two and three dimensions; in this paper we expand on 
these results both in depth and breadth for the case of 
steady state fiow. The results obtained reveal the rich 
and surprising nature of the collective behavior of the 
system. For certain values of the parameters, we observe 
Bagnold scaling in stable steady state fiow, with a con- 
stant density profile independent of depth. However, we 
also saw surprising examples of self-organization, includ- 
ing the fiow-induced crystallization of a disordered state 
into one with much lower dissipation. In this regime (sys- 
tems flowing on moderately smooth bottom surfaces) , we 
found reentrant disordering as well, and even oscillations 
between ordered and disordered states. The effects of 
bottom surfaces are thoroughly discussed in a separate 
work [^. In this paper, we concentrate on rough bottom 



surfaces for which the behavior is simpler. 

These simulations also allow us to investigate more 
subtle aspects of chute flow, such as hysteresis in the an- 
gle of repose and normal stress inequalities not accounted 
for by any conventional continuum theory. Additionally, 
we were able to look for surface and bulk instabilities to 
flow at the angle of repose. In particular, we find that 
although the Bagnold rheology of flow near the angle of 
repose is a bulk rheology, the fundamental instability in- 
ducing the flow in three dimensions appears to be an 
instability of the surface layers of the granular medium. 

Because granular materials are so common in nature, 
existing on many different length scales, there is a great 
amount of experimental data on a wide range of dynamic 
situations; shear flow and vibration experiments 
and studies of geological debris flows to name just 
a few [^0. There have been several moderately well- 
characterized experimental studies of granular flow down 
an inclined plane under laboratory conditions JlT|-p^. 
Yet for all the intense activity in this field over the years, 
the rheology of granular systems still remains a largely 
unsolved problem. 

There has been some work on continuum modeling of 
chute fiow; for a review of continuum based ideas see Sav- 
age and references therein. Other theoretical anal- 
yses (sometimes combined with case-specific simulation 
verification), specifically applied to chute flow geome- 
tries, attempt to calculate density and velocity profiles 
p8|-pT|, but a general consensus on the qualitative fea- 
tures of these profiles has yet to be reached. 

The state of the art of computer simulations of chute 
fiow is much less satisfactory because of the enormous 
equilibration times needed to set up steady fiow. In 
three dimensions, simulations have been performed for 
rather thin piles, which provides insight into only a small 
region of phase space | p^ , p^ , p3[ . Simulations of two- 
dimensional flows also report on small systems, and it 
is unclear whether these studies are carried out in the 
steady state regime or whether the data reported are 
transient |2^-p7[|. The basics behind granular simula- 
tions are available in Refs. ||2^ for 2D and for 3D. 

Our simulations attempt a systematic 3D study of 
chute flows. Unfortunately, we probe regions of phase 
space difficult to access in experiment. In a typical 3D 
experiment the flow is induced through a hopper-feeder 
mechanism which controls the flow rate of the system, 
but not the thickness of the flowing sample, which is 
chosen spontaneously by the system. Thus, much ex- 
perimental data is for flowing piles 10-15 particles high, 
whereas most of our simulations focus on moderate to 
thick piles, greater than 30 particles. Simulation results 
for systems smaller than 10 - 15 particles high do not 
show the same scaling as that for thicker systems 
Also, 3D experiments are usually carried out in narrow 
channels of order 10 particles wide or less, where side- 
wall effects may have a significant role in the observed 
behavior. Our simulations are periodic in this, the vor- 
ticity, direction, and we have yet to study wall effects. 



We suspect that most discrepancies that may exist be- 
tween different experimental and simulation studies are 
due to such differences in the detailed nature of the sys- 
tems studied. 

Because of the complexity of flowing granular systems, 
it is useful to first define the region of study. In order 
to determine the phase boundaries of fully developed, 
steady state flow, we have performed a series of simu- 
lations of inclined plane gravity driven flows in two and 
three dimensions in an attempt to define the region of 
phase space for which steady state fiows exist. A typical 
configuration snapshot in 3D, defining the computational 
geometry, is shown in Fig. 0. 




FIG. 1. Typical 3D snapshot for chute flow: N=24,000, 
with bottom surface dimensions 20x20 diameters shown by the 
black particles fixed to the bottom plate; tilt angle 6 = 24°, 
coefficient of restitution e = 0.88, and friction coefficient 
fj. = 0.50. Flow IS directed down the incline. 

In our simulations, initiation of flow is achieved by tilt- 
ing at a large angle (24—30°) to induce fiow. This proce- 
dure removes any configuration construction history ef- 
fects. We then reduce the inclination to a lower angle 
and allow the simulation to run until we observe a steady 
state flow regime (if one exists) . We define steady state as 
fiow wherein the energy input from gravity balances that 
dissipated from friction and collisions; so that the total 
kinetic energy of the system reaches a macroscopically 
constant value. In this case, the results are independent 
of sample history. 

In Fig. 0, we draw phase boundaries for both two and 



three dimensional flows as a function of the external con- 
trol parameters: tilt angle 6 and pile height H. This 
should be compared to a similar experimental determi- 
nation recently obtained by Pouliqucn |^^. The salient 
features are the existence in both 2D and 3D of three 
principal regions, corresponding to no flow, stable flow, 



and unstable flow. For a system of given thickness and 
fixed microscopic interaction parameters, these three re- 
gions are separated by two angles: Or, the angle of repose, 
and 9 max, the maximum stability angle, the largest angle 
for which stable flow is obtained, shown by a solid and 
dashed lines in Fig. ^, respectively. 




e 

(a) (b) 

FIG. 2. Phase behavior of granular particles in chute flow geometry, characterized by pile height H vs. tilt angle 9 for 
monodisperse systems in (a) 2D with fi = 0.50 and e = 0.92 (identified as Model L2 in Table and (b) 3D with fj. = 0.50 
and e = 0.88 (identified as Model L3). Both figures are for the spring dash-pot interaction model with rough bottom surface. 
Solid circles indicate region of steady state flow, open symbols correspond to no flow or unstable flow. In 3D we have identified 
hysteretic flow as x. 



For < Or, granular flow cannot be sustained. In 
the region Or < < Omax, we obtain steady state flow 
with packing fraction independent of depth. The region 
of constant packing fraction in the flowing material for 
steady state systems is accompanied by a smoothly vary- 
ing, non- linear velocity proflle. For > Omax, the devel- 
opment of a shear thinning layer at the bottom of the 
pile results in lift-off and unstable acceleration of the en- 
tire pile. The exact locations of these phase boundaries 
depend on the model parameters such as fi and e. For 
instance, in 2D, if e is reduced from 0.92 to 0.82, the max- 
imum angle of steady-state flow increases from « 23° to 
26°. Similarly, reducing fj, typically reduces the range of 
stable flow. 

It is well known that granular systems exhibit hys- 
teresis. Such behavior is usually attributed to system 
preparation and associated history effects . Although 
we observe three distinct regimes, the behavior close to 
the phase boundaries is sensitive to the procedure for the 
initiation of ffow. Indeed, we have observed hysteresis in 
our 3D simulations when approaching Or from either side, 
particularly for thinner piles. The hysteresis was signif- 
icantly reduced upon increasing pile height H. Crystal- 



lization of the 2D monodisperse pile upon the arrest of 
flow was primarily responsible for the large hysteresis ob- 
served in that case. 

Besides the phase diagram, our most important results 
concern the detailed structure and rheology of the steady- 
state flowing regime. In this regime, we do see a constant 
density proflle with height, as well as the Bagnold scal- 
ing of Eq. (|l]) . The amplitude of the strain rate goes to 
zero at the angle of repose; thus relations such as Eq.(|^) 
possess an additional strong angular dependence. 

We also analyzed the normal stresses in the flowing 
state, and found a number of results, most notably that 
the normal stress perpendicular to the free surface, azz, 
is approximately, but not exactly, equal to the normal 
stress parallel to the flow, axx- 

There are two fundamental puzzles in these results for 
the rheology of chute flow. The flrst, and smaller, puzzle 
is the appearance of an anomalous normal stress differ- 
ence azz — o'xx- We have been unable to define a simple, 
local, dimensionally consistent and rotationally invariant 
constitutive relation connecting a to 7^ that recovers this 
behavior. 

The second, and deeper puzzle, is the relationship be- 



tween the rheology and the Coulomb yield criterion. As 
the angle of repose is approached from above, the ampli- 
tude of the flow goes to zero, but the tensor structure of 
a remains approximately liquid-like, instead of recover- 
ing the large normal stress difference characteristic of the 
Coulomb yield criterion, which presumably applies to the 
static pile at the angle of repose. The one exception to 
this observation is the surface stress in three dimensions, 
where the normal stress differences do become large as 
the angle of repose is approached, suggesting that sur- 
face yield may control the failure of the static state. 

In the bulk, however, we are left with a transition to 
a static state that appears continuous in shear rate, but 
is apparently discontinuous in normal stress. We do not 
believe that an understanding of chute flow rheology is 
possible without resolving this seeming paradox. 

We present the simulation scheme in Section ||, detail- 
ing the inter-particle force laws. In Section^, we report 
our comprehensive simulation analysis, including the be- 
havior of the density and velocity profiles for our systems 
with varying interaction parameters. In Section [V 



we 



present a detailed discussion of stress analysis and rhe- 
ology of chute flow systems. In Section ^ we summarize 
our findings. 



II. SIMULATION METHODOLOGY 

We use the methods of molecular dynamics to perform 
2D and 3D simulations of granular particles. For this 
study we model N mono-disperse spheres of diameter d 
and mass m, supported on the xy-plane by a rough bed. 
The computational geometry of the present 3D system 
consists of a rectangular box with periodic boundary con- 
ditions in the x (flow)- and y (vorticity)-directions and 
constrained in the vertical z-direction by a flxed rough, 
bottom wall and a free top surface, as in Fig. |l|. Simu- 
lations in periodic cells attempt to study flow down in- 
finitely long and wide chutes, while using a finite number 
of particles. 

In 3D, the fixed bottom is constructed from a random 
conformation of spheres of the same diameter d as those 
in the bulk by taking a slice with areal fraction very close 
to random close packing (approximately 1.2 particle di- 
ameters thick), from a previously packed state. This sim- 
ulates an experimental procedure whereby glue is spread 
over the original smooth chute surface and particles are 
then sprinkled onto this surface to construct a rough floor 
approximately one particle layer thick. For 2D studies, 
the bottom wall is constructed from a regular array of 
spheres of diameter 2d and particle motion is restricted 
to the xz-plane. 

We employed a contact force model that captures the 
major features of granular interactions. In 2D, inter- 
actions between (projected) spheres are modeled using 
a linear spring model with velocity-dependent damping 
(the spring-dashpot interaction) and static friction. In 



3D the spring-dashpot model and static friction are also 
used, as well as Hertzian contact forces with static fric- 
tion. In the presentation of the results, we will specify 
which model is employed, and discuss the differences. 

The implementation of the contact forces, both the 
normal forces and the shear (friction) tangential forces, 
is essentially a reduced version of that employed by Wal- 
ton and Braun |22|, developed earlier by Cundall and 
Strack |28| . More recent versions of these models now ex- 
ist [pfl-Ba. We ignore hysteretic effects between loading 
or unloading normal contacts and we do not differentiate 
between frictional directions at the same contact point 
at different time-steps as does Walton [ p3| . 

Static friction is implemented by keeping track of the 
elastic shear displacement throughout the lifetime of a 
contact. For two contacting particles {«, j}, at positions 
{ri,rj}, with velocities {vi,Vj} and angular velocities 
{wijWj}, the force on particle i is computed as follows: 
the normal compression 5ij , relative normal velocity v„ . . , 
relative tangential velocity vj .j. are given by 



(2) 
(3) 

(4) 



where Ty = - rj, n,y = y^ij/nj, with r,y = |ry |, and 
Vy = Vj — Vj . The rate of change of the elastic tangen- 
tial displacement Ut , set to zero at the initiation of a 
contact, is given by 



du 



■tij 



dt 



(5) 



The second term in Eq.(||) arises from the rigid body 
rotation around the contact point and insures that u^.^. 
always lies in the local tangent plane of contact. Normal 
and tangential forces acting on particle i are given by 

Fn,j = f{Sij/d){knSijnij - 7„mcffV„^J, (6) 
Ft,, = f{d/d){-ktut^^ - 7tmeffVt,J, (7) 

where kn^t and jn,t are elastic and viscoelastic constants 
respectively, and nics = mirrij / {rrii + rrij) is the effective 
mass of spheres with masses rrii and rrij . The correspond- 
ing contact force on particle j is simply given by Newton's 
third law, i.e., Fji = — . For spheres of equal mass m, 
as is the case here, rrics = m/2; f{x) = 1 for the lin- 
ear spring-dashpot model, denoted henceforth as Model 
L, or f{x) = y/x for Hertzian contacts with viscoelastic 
damping between spheres, denoted as Model H. 

Our results are given in non-dimensional quantities 
by defining the following normalization parameters: dis- 
tances, times, velocities, forces, elastic constants, and 
stresses are respectively measured in units of d, to = 
■y/ d/g, Vo = Vgd, Fo = mg, kg = mg/d, and Uo — 
mg/d^. For a realistic simulation of glass spheres with 
diameter d = lOO/im, the appropriate elastic constant 
j^giass _ q(^iqW^ necessitates a very small time-step for 



accurate simulation, prohibiting any large-scale study. 
In our simulations, we typically use a value for fc„ = 
O(IO^) which we believe captures the general behavior 
of intermediate-to-high-A: systems, thus offering a rea- 
sonable representation of realistic granular materials (we 
discuss this aspect further in Section IIIB). A complete 
list of model parameters used in our standard simulation 
set, which consists of 2D and 3D versions of Model L (L2 
and L3), and a 3D version of Model H (H3), are given in 
Table |. 

In a gravitational field g, the translational and rota- 
tional accelerations of particles are determined by New- 
ton's second law, in terms of the total forces and torques 
on each particle i: 



(virial) and kinetic terms, 



j 



(8) 
(9) 



The amount of energy lost in collisions is characterized 
by the inelasticity through the value of the coefficient of 
restitution. For Model L, there are separate coefficients, 
e„ and et , for the normal and tangential directions which 
are related to their respective damping coefficients j„^t 
and spring constants kn^t- 



where the collision time tcoi is given by. 



(10) 



(11) 



The value of the spring constant should be large enough 
to avoid grain interpenetration, yet not so large as to re- 
quire an unreasonably small simulation time step 5t, since 
an accurate simulation typically requires 5t ^ tcoi/50. 
For Model H, the effective coefficients of restitution de- 
pend on the initial velocity of the particles. 

The static yield criterion, characterized by a local par- 
ticle friction coefficient /i [Q, is modeled by truncating 
the magnitude of u^.^ as necessary to satisfy IF^.^. | < 
|/xF„;^- 1. Thus the contact surfaces are treated as "stuck" 
while Ft . . < fiFm , and as "slipping" while the yield cri- 
terion is satisfied. This "proportional loading" approxi- 
mation is a simplification of the much more compli- 
cated and hysteretic behavior of real contacts To 
test the robustness of the proportional loading assump- 
tion, we also carried out simulations with Model L in 
which U( is not truncated but the local yield criterion 
Ft < fiFm is implemented. Note that we do not believe 
this to be a physically reasonable choice. Results for the 
two cases are similar, although the average kinetic energy 
is somewhat smaller (by approximately 18% for Model 
L2) when ut is truncated compared to those simulations 
when Ut is unbounded. 

The components of the stress tensor aap within a given 
sampling volume V are computed as the sum over all par- 
ticles i within that sampling volume of the contact stress 
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V 



ij 



(12) 



where F^ = F^. _ +Ff.. , and v is the time-averaged veloc- 
ity of the particles within the sampling volume V. The 
time-averaged velocity must be subtracted since the ki- 
netic portion of the stress tensor is entirely due to fiuc- 
tuations in the velocity field. 

For Hertzian contacts the ratio kt/kn depends on 
the Poisson ratio of the material, and is about 2/3 for 
most materials. For ease in our simulations, we use a 
value fct/fc„=2/7, which makes the period of normal and 
shear contact oscillations equal to each other for Model 
L ll^. However, the contact dynamics are not very sen- 
sitive to the precise value of this ratio. We have per- 
formed simulations with different values of kt / kn to test 
how this ratio may affect our results; different values of 
this ratio yield nearly identical results. The only differ- 
ence we observe is a slight increase in the total, averaged 
kinetic energy (KE) of the system when kt/kn > 2/7, 
and a decrease for kt/kn < 2/7. For example, when we 
set kt/kn — 2/3 instead of 2/7, the total averaged KE 
increases by about 10%, whereas all other macroscopic 
quantities measured in the simulations, such as density 
and stress, remain essentially unchanged. 

Similarly, although all results reported here are for 
7t/7n = {i.e. no rotational velocity damping term) we 
have also carried out simulations to measure the effect 
of introducing rotational damping, jt/jn > 0. When we 
set 7t = 7„, we observe a slight decrease, of about 8%, in 
the total averaged KE, compared with those simulations 
that have "ft = 0. Making jt/^n non-zero quickens the 
approach to the steady state by draining out more energy. 
However, all other quantities are, again, unchanged. We 
discuss reasons why we observe minimal changes with 
these interaction parameters in Sec. 



IHB 



Typical values for the friction coefficient fi range be- 
tween, 0.4 and 0.6 for many materials. We chose fj, = 0.50 
for most of our simulations, though variations in fi will be 
discussed in Sec. [II B. Similarly, the value of e is chosen 



to reflect the properties of a realistic granular particle. 
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TABLE I. Parameter values used in our standard simula- 
tion set for the 2D and 3D linear spring models (Models L2 
and L3, f{x) = 1), and the 3D Hertzian model (Model H3, 
f{x) = \fx). For Model H3, t is velocity dependent. 



The equations of motion for the translational and ro- 
tational degrees of freedom are integrated with either 
a third-order Gear predictor-corrector or velocity- Verlet 
scheme Q with a time-step St — 10^** for — 2 x 10^. 
All data was taken after the system had reached the 
steady state. To reach the steady state, simulations were 
required to run for 1 — 2 x lO"^ St when starting from a 
non-flowing state for N < 10, 000, and the largest system 
in 2D {H = 200) requhed a run time of 2 - 5 x lO^St. 
On a 500 MHz DEC Alpha processor, our code requires 
about 5 days to simulate 10 million timesteps of a 3D 
8000-particle granular system. We have also created a 
parallel version of the 3D code using the standardized 
MPI message-passing library. The parallel code parti- 
tions the simulation domain into small 3D sub-blocks us- 
ing the methods described in Even on a cluster 
computer with relatively low interprocessor communica- 
tion bandwidth, the code runs at high parallel efflciencies 
so long as we simulate 1000 or so particles per processor. 
For example, on 8 processors of our Alpha/Myrinet clus- 
ter, we can simulate 15 million timesteps/day of the same 
8000-particle system. 

For the imposition of chute flows with varying tilt an- 
gles, we rotate the gravity vector g in the xz-plane by the 
tilt angle 9 away from the — z direction — flow is from 
left to right in this sense. This means that the z-axis 
is always normal to the free surface. In 3D the area of 
the bottom is j4 = L^Ly where and Ly are the di- 
mensions of the simulation cell in the x and y directions 
respectively. For the 3D simulations, we define a mea- 
sure of the height of the pile by defining H = NcP/A as 
the pile height if it were sitting on a level plane at rest 
in a simple cubic lattice. For example, for iV = 8000 
and Lx — 2Qd and Ly — IQd, H — 40 (although due 
to the precise configuration, the actual measured height 
« 37). This is a useful definition for comparing between 
different system sizes. We study a range of system sizes, 
1000 < N < 20000. For the largest system, H = 100. 
The influence of other wall dimensions L^, Ly was also 
studied. For the 2D runs, the cc-dimension of the periodic 
side is flxed at lOOd {i.e. 50 large particles long) and the 
pile height 2<H< 200, i.e. N = 200 - 20000. 

In 2D, the initial state was constructed by building a 
triangular lattice of particles. The tilt angle was then 
increased until flow occurred. The initial flow occurred 
only for 9 ^ 23°. This minimum value to induce flow de- 
pends on the size and spacing of bottom plate particles. 
The initial failure occurred mostly at the bottom of the 
pile, followed by movement of a dilation front toward the 
top of the pile as shown in Fig. ||. Once this initial steady 
state was achieved, the angle 9 was adjusted to its de- 
sired value, and the system equilibrated to its final steady 



state. In 3D, we started the system from a randomly di- 
luted simple cubic lattice. The angle was then increased 
to a large angle 9 — 30° to induce disorder and settling of 
particles. The angle 9 was then decreased to the desired 
value and flow allowed to continue until a steady state 
was reached, before measurements were taken. 

In 3D, to test for hysteresis near 9^, 9 was reduced to 
below 9r until the system settled down into a disordered 
state and stopped flowing. 9 was subsequently increased 
to 9^'°"' and the system began flowing. This angle of 
flow initiation was sometimes different from the angle of 
cessation of flow 9f°P when taking a flowing state and 
then lowering 9 down to 9f"P to stop the flow. However, 
this small hysteretic behavior, in 3D, only occurs for thin 
piles at low angles. 

In 2D the equivalent phase diagram can only be con- 
structed by taking a flowing state at angle 9 and then 
lowering to 9f°P. Once the 2D state stops flowing the 
system spontaneously crystallizes into a polycrystalline 
ordered state. To induce flow from this ordered state 
requires increasing 9 to a much higher angle than 9f°P. 

III. RESULTS: VELOCITY AND DENSITY 
PROFILES 

A. Kinematics of steady state systems 

We focus our main attention on the regime of steady 
state flow for moderate to deep piles, for which 9r is in- 
dependent of depth. In Fig. Q we plot the density and 
velocity (in the direction of flow) z-proflles over a range 
of inclination angles 9, for a series of simulations in 2D 
and 3D. Figure |(a) is for a 2D system (Model L2, cf. Ta- 
ble |) of N=10,000 particles, corresponding to H = 100. 
In Fig. ^(b), the equivalent 3D model (Model L3 with 
N = 8,000, — 40), denoted by the open symbols, is 
compared to the 3D Hertzian model (Model H3). The 
tilt angle was varied between 18° — 30° in all cases. In 
2D the system becomes unstable to an accelerating flow 
above 23° and in 3D the unstable flow regime is observed 
above 26°. 

In both 2D and 3D, the packing fraction remains con- 
stant over almost 40 layers in 3D, and 100 layers in 2D. 
For all steady state systems, as the tilt angle is increased, 
the value of the bulk packing fraction decreases. This 
decrease accompanies a growing dilated region (of lower 
packing fraction) near the free surface at the top. All 
the velocity proflles are concave, and velocities increase 
in value with increasing tilt angles. Consequently, the 
total kinetic energy of the system rises with increasing 
angle [||. 





FIG. 3. Time sequence of a typical configuration in 2D following an instantaneous change in the inclination angle 6 from 0° 
to 24°. Results are for N = 10000, fi = 0.50, and e = 0.82, and for times t = a) 100, b) 400, c) 600, d) 6000. Flow is left to 
right. As the flow progresses, the dilational front propagates upwards through the system, destroying the initial ordered array; 
the pile consequently "fluffs" up. 




FIG. 4. Packing fraction tf) and velocity profiles, as a function of distance from bottom z, for (a) 2D spring-dash-pot model 
(Model L2), with H = 100, at tilt angles of 9=18, 19, 20, 21, 22, 23 degrees, (b) 3D, H = 40 systems at 9 =20, 22, 24, 26 
degrees, with Model L3 (open symbols) and Model H3 (solid lines). 



We monitor vertical mixing of the bulk by measuring 
the bulk-averaged, mean-square displacement of particles 
over time. Figure |5| shows the mean-square displacement 
of particles normal to the surface < > as a function 
of simulation time for Model L3, over a range of tilt an- 
gles. The linear relationship demonstrates well-defined 
diffusive motion in the z-direction, suggesting thorough 
mixing in the system. Similar results are observed in 2D. 
At long times < > will reach a constant due to the 
finite height of the pile. 
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FIG. 5. z-component of the mean square displacement for 
3 angles, 9 = 21°, 22°, and 23°, for Model L3, with H = 40. 

By observing a sequence of snapshots (not shown here) 
of tracer particles at various heights in the bulk, we also 
find that diffusion is somewhat faster nearer the bottom 



of the pile. This is indicative of the fact that fluctuations 
in the particle velocities are greater closer to the bottom 
wall. Fig. H depicts the diagonal components of the ki- 
netic part of the stress tensor, p < {5v°'Y >i where p is 
the mass density and Jw" = v°- — v°', at three different 
angles for Model L3. Indeed we do find that the velocity 
fluctuations (frequently termed "granular temperature" 
in the literature of dilute granular flows) are greatest at 
the bottom of the pile (away from the actual plate) and 
decrease with height until the values appear to level off 
at the top free surface. 
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FIG. 6. Profiles of the kinetic portion of the diagonal el- 
ements of the stress tensor p < (5«")^ >, normalized by 
their maximum value along the curve, for Model L3 inclined 
at 21° ( ), 23° (— ), and 25° {—). 



This behavior partially illustrates how the pile is able 
to maintain a constant density profile, even though the 
stresses increase towards the bottom, and the flowing pile 
has a finite compressibility, as evidenced by the chang- 
ing density as a function of tilt angle 9. Particles deeper 
into the pile experience increasing compaction forces due 
to the load of the particles above, yet a constant den- 
sity is maintained through the increased particle velocity 
fluctuations. 

The data sets shown in Fig. ^ are for one system size 



only. In Fig. ^ density and velocity profiles for systems 
of varying heights are compared. The densities measured 
deep in the pile, as well as the density and strain rate 
profiles near the surface, are independent of the overall 
height of the pile. This suggests that the rheology of the 
system is local in this regime; i. e., that constitutive re- 
lations locally rel ate st ress and strain rate. For reasons 
alluded to in Sec. [VD, we have been unable to identify 



these constitutive relations. 
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FIG. 7. Density and veloctty profiles for (a) 2D systems (Model L2) for 6 = 20° with sizes H = 200, 100, 50, 25 and (b) 
3D systems (Model L3) for 6 = 24° with sizes H = 100, 60, 50, 60 . 
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FIG. 8. Density and velocity profiles for thin systems 
(Model L2), with H = 10, at 20° and 23°. These profiles 
are very different from the thicker piles. 



We note that the behavior observed in Figs. ^ and |^ 
is true only for H > 20. For smaller piles, such as those 
commonly studied experimentally jl^ , the behavior is 
very different, as seen in Fig. |^, where a linear velocity 
profile is evident. Recent experimental studies of thin 
piles in inclined plane geometries also report linear ve- 
locity profiles [ pl| . 

B. Dependence on Interaction Parameters 

In this subsection, we investigate the sensitivity of 
these results to the particle interaction parameters. We 
independently vary the internal coefficient of friction fi, 
the coefficient of restitution e, and the value of the spring 
constant fc„. We observe that while the density of the 
bulk material does not depend sensitively on these inter- 
action parameters, the velocity profiles do. 

Figure ^ shows the sensitivity to the friction coefficient 
fi by depicting density, velocity, and strain rate profiles 
for: (a) Model L2 with = 50 and 6* = 20°, where 
fi = 0.15, 0.25, 0.50, 1.0, and (b) Model L3 with H = 40 



and e = 22°, for values of ^ = 0.15, 0.25, 0.5, 1.0. The 
data suggest that there is minimal change in the bulk 
density over this range in ^. In the bottom panels of 



Fig. ^ the shear rate scaled by 
the various values of ^. 
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FIG. 9. Density, velocity, and strain rate profiles for different values of the particle friction coefficient for: (a) Model L2 at 
= 20°, = 50 /or = 0.15, 0.25, 0.50, 1.0, and (b) Model L3 at 6 ^ 22° , H = AO for fj, = 0.15, 0.25, 0.50, 1.0. 



Similarly, Fig. shows the profiles for the same sys- 
tems as described in Fig. ||, but with a fixed fi — 0.5 and 
varying coefficients of restitution e, i.e., varying the in- 
elasticity of the system. Again we see that variations in 
e have little effect on the flow behavior of these systems. 



particularly in 3D, provided that the system is able to 
reach steady state. [For low ^ (« 0.10) and high e (« 
0.96 for 2d and 0.98 for 3d), the systems become unsta- 
ble] 




FIG. 10. Density and velocity profiles for different values of e for (a) Model L2 at 9 = 20° , H = 50 for e = 0.72, 0.82, 0.92, 0.96, 
and (b) Model L3 at 9 = 22°, H = 40 for e = 0.58, 0.78, 0.88, 0.98. 



Another microscopic parameter we have investigated 
is the effective hardness of the particle, determined by 
the value of the spring constant /c„. We vary fc„ and 
keep e constant by adjusting the value of 7„. Simulations 
investigating this parameter can be time-consuming: in- 
creasing kn by a factor of 100 requires a reduction in 
the time-step by a factor of 10. Fortunately, as Fig. [Tl| , 
(measured for Model L2 with 9 = 20°, H = 50) indicates, 
the effect of variations in kn is minimal, provided kn is 
sufficiently large. 
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C. Dependence on Tilt Angle 

Judging by the insensitivity of the macroscopic quan- 
tities to the various interaction parameters for Model L 
(as shown), as well as Model H, we see that to a good ap- 
proximation, effects due to material properties and sys- 
tem size can be neglected in the steady state regime. 
As shown in Fig. |l^, the packing densities vary approxi- 
mately linearly with 6 and approach the maximum values 
(j)fg^ = 0.815(5) and ^g^^'^ = 0.590(5) at 6'^,2r> ~ 17.8° 
and 6r,3D ~ 19.4° for 2D and 3D, respectively. 

It is interesting to note that the asymptotic packing 
fractions and <j)^^'^' are close to the values one would 
obtain assuming the flow was the densest possible flow 
of lines (in 2D) or planes (3D) of close-packed particles 
parallel to the top surface. For the 2D case, the packing 
corresponds to a square lattice with a packing fraction of 
7r/4 ~ .79. For the 3D case, the sliding planes would be 
square lattices, stacked to form triangular lattices in the 
y — z plane. This arrangement has a packing fraction of 



IV. RESULTS: STRESS ANALYSIS 



60 



A. Cauchy Equations 



FIG. 11. Density and velocity profiles for different values 
of the spring constant for Model L2 at — 20°, H = 50. 
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FIG. 12. Tilt dependence of the packing fraction in the re- 
gion of constant packing fraction, for Models L2 (solid cir- 
cles), L3 (solid squares), and H3 (open diamonds). The 
dashed lines denote the linear dependence on tilt angle near 
the angle of repose. 



The stress tensor is symmetric: ct^ = Uji, with 
D{D -\- l)/2 independent components in Z)-dimensional 
space. The Cauchy (force-balance) condition provides 
only D equations, leaving the solution underdetermined. 
Thus, an additional D{D — l)/2 constitutive relations 
are needed to close the equations and to solve for the 
transmission of stress in a granular system. 

In 2D, the steady state Cauchy equations are 



dz 
daxz 
dz 



pg cos 9, 



pg sm ( 



For a given tilt angle, these give: 



,{z) — g cost 



dzp{z) 



(Txz{z) = azz{z)ta,ii9, 



(13) 
(14) 

(15) 
(16) 



where p is the number density of spheres {4> — npd^ /2D 
for dimensionality D = 2 and 3). If, as in our case, the 
density p is constant. 



cTzziz) = gpcose{h - z), 
<7xz{z) = gpsm9{h - z), 



(17) 
(18) 



where h is the effective height of the flowi ng p ile, which 

appears as a constant of integration in Eq. ( |l7| ) . 

not be determined from these considerations; as we lack 



a constitutive relation that would determine it. Never- 
theless, important features of the behavior of the stress 
tensor can be obtained by Mohr-Coulomb analysis 



if the flowing pile behaves like a fluid, <Tx 
consequently sin2iy9 = tan0. 



and 



B. Mohr-Coulomb Analysis 

The Mohr circle, shown in Fig. [l^a, is a geometri- 
cal construction that enables visualization of rotational 
transformations of the stress tensor. The circle is drawn 
in the cr — r plane, such that the points A {uzz, Cxz) and 
B (axx, —Cxz) form a diameter of the circle, centered at 
point O. Coordinates of the points on the circle represent 
the normal (a) and shear (r) components of the stress 
tensor associated with all possible shear planes. Upon a 
rotation of the coordinate system, i.e., the plane in which 
shear is specified, by an angle tp, the representative points 
rotate by an angle 2ip around the circle. 

At a given tilt angle, azz and axz are determined by 
the Cauchy equations, which fixes the location of point A 
{<^xz/<^zz = tan^). However, axx, and thus the location 
of point B, is undetermined by the Cauchy equations, 
and depends on the rheology. 

The "stress angle" 2ip = CO A, formed by the stress 
point A, origin of the Mohr Circle O, and point C, whose 
tangent passes through the origin of the (ct, r) plane, can 
be used as a surrogate for any quantity that completes 
the description of the x ~ z stress state, since it uniquely 
identifies the two-dimensional stress state of the flowing 
pile (for 9 > Or) by fixing the value of axx- For a pile 
with a uniform Coulomb yield criterion that is at incip- 
ient yield everywhere (lYE) when 6 = 9r, the points C 
and A coincide, and therefore tp = 0. On the other hand. 



C. Stress Tensor Near the Surface 

In all cases, the behavior of 2(p as a function of depth 
can be fitted to an empirical form that starts at a "sur- 
face" value at the effective height h and approaches a 
"bulk" value exponentially (see Fig. p^): 

Figure ^ depicts the values for the fitting parameters 
2lp^^'^^ and 2lp^^^^ as a function of tilt angle for the three 
main models studied in this paper. The following obser- 
vations can be made: 

(i) 2ip, and consequently all the ratios of stress ten- 
sor components, becomes independent of depth below a 
transitional surface layer about 5d to 8d in thickness. 

(ii) In 2D (Model L2, see Fig. |I|a), as 6 is lowered 
to 6r, the stress state at the surface moves even farther 
from lYE compared to the bulk. Independent observa- 
tions confirm that the top surface does not play any dis- 
cernible role in the arrest and start of flow; this primarily 
occurs near the bottom surface. 



(iii) However, for both models in 3D (Fig. 143), as 9 is 
lowered to 9r, the surface layer does approach incipient 
yield {2ip = 0) while the bulk remains far from it. It ap- 
pears that the stabilization of the surface layer aX 9 = Or 
is responsible for the arrest and subsequent restart of flow 
in the entire system, accompanied by a near-elimination 
of flow hysteresis. 




(a) (b) 

FIG. 13. (a) The Mohr circle is a graphical tool that is used to determine transformations of a rank 2 tensor (such as stress) 
under rotation. The stress components for a given coordinate system are represented by points A and B, which form a diameter 
of the circle. The transformed stress components upon a rotation of the coordinate system by angle ip can be found by a rotation 
of these points by 2tp around the circle. The point C, which has a tangent that passes through the origin, corresponds to the 
orientation of a shear plane (at an angle tp to the x—axis) with the largest ratio of shear to normal stress, (b) The stress angle 
2ip [CO A in (a)] as a function of height for 9 = 20° (Q ), 22° (a ), 24° (O ), and 26° (A ). The results are for Model H3 with 
H=40. The lines are fits that decay exponentially from 2(/5'^"'^ at the effective height h [cf. Eg. ^P\)] to 2</p'°""' in the bulk, with 
a typical decay length of 1.3 to 2.2d, indicating a surface layer about 5d to 8d thick. 
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FIG. 14. The stress angle at the surface, 2<^°"'^^ (circles), and in the bulk, 2;^''""^ (squares), for (a) Model L2 (open symbols 
connected by dotted lines), and (b) Model L3 (open symbols connected by dotted lines) and H3 (solid symbols connected by dashed 
lines). For Model L2, the rheology at the surface (2 = h) near 6r is even farther away from the lYE condition compared to the 
bulk. However, both 3D models observe near-IYE conditions at the surface near 9,., suggesting that the arrest of flow may be 
initiated by the surface rather than the bulk. The solid lines depict behavior expected without a normal stress anomaly, when 
sm2(fi = tan 8. 



(iv) The bulk has nearly identical normal stresses 
axx and azz, which would have corresponded to 2(p = 
arcsin(tan0), depicted by the solid lines in Fig. In 
other words, the normal stress anomalies discussed in 
Sec. p^VD are quite small compared to what one would 
have attributed to a plastic material at incipient yield. 

(v) The transitional surface layer is not directly re- 
lated to the dilated layer; the former is much thicker near 
6 = Or and penetrates well into the region of constant 
density, as can be seen by comparing Figs. ^ and |l3|b. 
In fact, upon approaching 9^, the width of the surface 
rheological layer 6 increases slightly whereas the width 
of the dilated layer decreases substantially. 
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FIG. 15. Profiling the components of the stress tensor at 
9 = 22° m (a) Model L2 for H = 100, (b) Model L3 for 
H = 40. 



D. Bulk Rheology 



Having identified the behavior associated with the free 
surface at the top, we can now investigate the stress 
tensor below this surface layer. For tilt angles suffi- 
ciently above Or, where the granular medium behaves 
roughly like a fluid, one might expect the normal stresses 
{axx, o'yy, azz in 3D, axx and azz in 2D) to be equal. 
In 3D, we find that ayy is smaller than the xx— and 
zz— components by 15 — 20%, suggesting that consolida- 
tion and compaction normal to the shear plane is poorer. 
The normal stresses and the driving shear stress axz , for 
the 2D and 3D linear spring model are shown in Fig. |l^. 
The components {axy, ayz) are not shown, since they van- 
ish due to the symmetries in the geometry; they are in- 
deed measured to be zero within the error bars associated 
with the sample size and averaging time. 

We observe that for both the 2D and 3D systems 
(and for both linear spring and Hertz models), although 
a XX ~ azz , there are small but systematic deviations from 
perfect equality that become independent of depth in the 
bulk. Let us define a "normal stress anomaly" x as. 



(20) 



This is simply an alternate parameterization of the stress 
angle 2ip defined earlier, introduced as a convenience to 
emphasize the small deviations around fluid-like behav- 
ior, for which % = 0. Therefore, x is also independent of 
height z except near the top and bottom surfaces. With 
this in mind, we plot the bulk value of x vs. in Fig. |l^, 
noting a strong angle dependence, in which x is neither 
monotonic in 6 nor of a specific sign. We have evaluated 



a class of homogeneous, polynomial, rotationally invari- 
ant constitutive stress-strain rate relations, but have not 
been able to satisfactorily describe these rather peculiar 
normal stress anomalies. 
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FIG. 16. Dependence of the normal stress anomaly x on 
tilt angle 6 for Models L2 and L3 (closed symbols), and Model 
H3 (open symbols). 

The fact that the stress varies linearly with depth 
and our earlier observations of constant density suggests 
that the analysis relevant to our systems is that due to 
Bagnold Bagnold's collisional-momentum transfer 
analysis for granular systems works under the assump- 



tion of a constant density profile, resulting in stress pro- 
files that vary linearly with depth. The essence of Bag- 
nold's theory is a constitutive equation whereby the shear 
stress Uxz is proportional to the square of the strain rate 
7^ = {dvx{z) / dzf' , where Vx{z) is the velocity in the 
direction of flow at height z: 



(21) 



Combined with Eq.(|18|), and the no-slip boundary con- 
dition at z = 0, this results in a velocity profile of the 
form, 



pg sm ( 



1 - 



h- 



3/2' 



(22) 



From Fig. |r^, we observe that for the bulk of the flow, 
the relationship axz oc 7^ holds to a good approximation 
below the first 5-8 layers, and away from the bottom wall, 
for the 2D and 3D systems. We have fitted the "Bagnold" 
scaling with the dotted lines, the solid lines and symbols 
represent the simulation data |^^. The tilt dependence 
of the overall amplitude of the strain rate, ^Bag, is shown 
in Fig. |l^. In 3D, ^Bag continuously approaches zero at 
the angle of repose, whereas in 2D, there is a jump in 
this amplitude, consistent with the overall hysteretic be- 
havior. 




FIG. 17. Rheology curves of chute flow systems-shear strain vs. shear stress; a) 2D, (Model L2, H = lOOj and b) 3D 
(H = AQ): Model L3 (symbols) and Model H3 (solid lines). Bagnold scaling fits m the bulk are shown by dashed lines for Model 
L and dotted lines for Model H. 
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FIG. 18. Strain rate amplitude Asag = j/y'axz associated 
with the bulk Bagnold rheology for the three model systems. 
Whereas 3D amplitudes extrapolate to zero at 6r, there is a 
finite jump associated with the 2D amplitude at Or . 

Another way to test this scahng is by plotting the av- 
erage velocity < v"^ >i/2 as a function of H (which is 
proportional to h). The scaling in Fig. ^ shows that 



< >i/2(x H", where a = 1.52 ± 0.05. This result also 
agrees well with experiment If we rescale the data 
from Fig. ^, we find good agreement apart from the re- 
gion near the top surface where the density is no longer 
constant. This suggests that Bagnold's theory may pro- 
vide an approximate description of the bulk motion of 
our systems. In fact, Bagnold scaling is a generic dimen- 
sional result for the situation where the time scale of the 
system is only set by the inverse of the shear rate, as is 
the case here [^. 

Bagnold's original stress-strain rate relationship arises 
from a momentum transfer mechanism that is based on 
binary collisions. From the simulation data, we find that 
the dominant term in the stress is due to lasting contacts 
between particles, and the ballistic (kinetic) contribution 
to the stress is significantly smaller (about 1% of the total 
value). Thus, the success of the Bagnold scaling is based 
on the dimensional structure of the problem, rather than 
on the particular momentum transfer mechanism that he 
identified. 
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FIG. 19. Scaling of velocity m the direction of flow < v'^(z) >^/^ with system height H in; (a) 2D at 20° with e = 0.92, for 
two different values of /i, (b) model H3 at 24°. The slope of the lines indicate that < v"^ >^''^oc H" , with a — 1.52 ± 0.05. 



Another method to test the nature of collisions is to 
compute the average coordination number Zc func- 
tion of inclination angle. This data, for the 2D and 3D 
linear-spring systems, is shown in Fig. In a system 
dominated by binary collisions, one would expect Zc <C 1; 
this is clearly not the case for our system. The observed 
behavior is an increasing Zc as 6 approaches the angle of 
repose from above. Normalized this way for a static 2D 
triangular lattice with no free particles, the value would 
be 3. Similarly, for 3D static packings one might expect 
a value between 4 — 6 [H . 



butions to the kinetic term of the stress tensor do not 
play a significant role in determining the macroscopic 
quantities measured. It might then be argued that for a 
densely packed pile of stiff objects in motion, the overall 
time evolution of the system in the configurational phase 
space is primarily constrained and controlled by aspects 
of geometrical packing, rather than the specific form of 
the stiff force laws between particles or dissipation func- 
tions. This might be why the system is so insensitive to 
vari ations in the interaction parameters, as described in 
Sec.lmBl 



Because of these observations, we reason that contri- 
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FIG. 20. Averaged instantaneous coordination number 
Zc,as a function of tilt angle for Model L2 with H = 100, 
and Model L3 with = 40. 



V. CONCLUSIONS 

We have concentrated on the steady-state nature of 
chute flows, specifying first the region in phase space in 
which such flows can be observed, and second the struc- 
ture and rheology of these flows. 

A region of constant packing fraction is a generic fea- 
ture in our 2D and (two) 3D models, with only a small 
dilated layer at the free surface. Analysis of the ve- 
locity proflles has revealed that to a good approxima- 
tion, Bagnold scaling holds: the Bagnold velocity profile, 
Vx (x H^-^, and rheology, cr cx 7^, is reasonably verified 
away from the surface. This is in contrast with earlier 
simulations on chute flows, which indicated linear veloc- 
ity profiles. We argue that although this latter may be 
the case for small systems, such as flowing layers less 
than 20 particles high, steady flows of moderately thick 
systems are well-approximated by Bagnold scaling. 

Although the regime of Bagnold-like flow appears to 
dominate the system, we have found that deviations from 
this simple theory exist. The normal stress anomaly re- 
mains a mystery, and our flts to the stress-strain rate 
curves apply only away from the top and bottom sur- 
faces. We have also found that the transmission of stress 
in such dense flows is dominated by contacts, as opposed 
to binary collisions in Bagnold's analysis of dilute flows. 

Finally, we observe that the normal stresses in bulk 
flows do not approach a Coulomb yield criterion struc- 
ture at the angle of repose, despite the continuous dis- 
appearance of the shear rate at this threshold. The fact 
that Coulomb yield is approached at the surface for 3D 
flows hints at a special role for surface failure in this case. 

Our simulation code, both in its simple and parallelised 
versions, enables us to study large systems for very long 
time scales, and we continue to investigate some of the 
outstanding issues in this area. We will report elsewhere 



the differences between rough and smooth bottom sur- 
faces We will also go on to study 3D planar Couette 
flows, extending Ref. |^^, and will be reporting on this 
in the future. 
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